index

Author

Liz Loutrage

1. Data

2. Isotopic niches

Ellipses

  • ellipses at 40%
  • Δ \(\delta\)13C = 2.36‰
  • Δ \(\delta\)15N = 5.94‰
Code 
niche_plot_community <- isotope_data%>%
  mutate(species= gsub("_"," ", species))%>%
    mutate(species=recode(species, "Cyclothone"="Cyclothone spp."))

# colors selection 
nichecol_krill <- c("#E4A33A", "#F67451", "#D664BE", "#3DA5D9", "#94b3ae",
                    "#18206F", "#FD151B", "#049A8F", "#072AC8", "black", "purple", 
                    "#d193f7", "#d8c2ab", "#678FCB", "#A63A49", "#00547A", "#6B54A0")

# size of the ellipse 
p.ell<-0.40
 
# plot
ggplot(data = niche_plot_community, 
                         aes(x = d13c, 
                             y = d15n)) + 
  geom_point(aes(color = species, shape= taxon), size = 2) +
  scale_color_manual(values=nichecol_krill)+  
  scale_fill_manual(values=nichecol_krill)+
  scale_shape_manual(values= c(19, 3))+
  scale_x_continuous(expression({delta}^13*C~'\u2030'), limits = c(-21, -18.5)) +
  scale_y_continuous(expression({delta}^15*N~'\u2030'), limits = c(7, 14))+
  stat_ellipse(aes(group = species, fill = species, color = species), 
               alpha = 0.2, level = p.ell,linewidth = 0.7, type = "norm", geom = "polygon")+
  theme_light()+
  theme(legend.text = element_text(size=15),
        legend.title = element_text(size=15),
        axis.title = element_text(size=15),
        axis.text = element_text(size=15))+
  labs(shape="Taxon", col= "Species", fill="Species")

Code 
#ggsave("niches_community.png", path = "figures", dpi = 700, height = 10, width = 12)

Niche areas

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 14
   Unobserved stochastic nodes: 3
   Total graph size: 29

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 120
   Unobserved stochastic nodes: 3
   Total graph size: 135

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 57
   Unobserved stochastic nodes: 3
   Total graph size: 72

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 20
   Unobserved stochastic nodes: 3
   Total graph size: 35

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 97
   Unobserved stochastic nodes: 3
   Total graph size: 112

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 20
   Unobserved stochastic nodes: 3
   Total graph size: 35

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 26
   Unobserved stochastic nodes: 3
   Total graph size: 41

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 20
   Unobserved stochastic nodes: 3
   Total graph size: 35

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 41
   Unobserved stochastic nodes: 3
   Total graph size: 56

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 20
   Unobserved stochastic nodes: 3
   Total graph size: 35

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 43
   Unobserved stochastic nodes: 3
   Total graph size: 58

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 14
   Unobserved stochastic nodes: 3
   Total graph size: 29

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 60
   Unobserved stochastic nodes: 3
   Total graph size: 75

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 12
   Unobserved stochastic nodes: 3
   Total graph size: 27

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 20
   Unobserved stochastic nodes: 3
   Total graph size: 35

Initializing model

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 20
   Unobserved stochastic nodes: 3
   Total graph size: 35

Initializing model

Overlaps

Code 
# Prepare data 
mat_overlap_data <- isotope_data%>%
  filter(taxon=="Fish")%>%
  mutate(species= gsub("_"," ", species))%>%
  mutate(species=recode(species, "Cyclothone"="Cyclothone spp."))

# Arrange species by their d15n values 
mat_overlap <- mat_overlap_data%>%
  group_by(species)%>%
  mutate(mean_d15n=mean(d15n))%>%
  arrange(mean_d15n)%>%
  as.data.frame()

# ellipse at 40%
test.elp_community <- rKIN::estEllipse(data= mat_overlap,  x="d13c", y="d15n", group="species", levels=40, smallSamp = TRUE)
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
Code 
# Extract the area of each polygon
elp.area_community <- rKIN::getArea(test.elp_community)

# determine polygon overlap for all polygons
elp.olp_community <- rKIN::calcOverlap(test.elp_community)%>%
  #renames rownames
  mutate(OverlapID= gsub("_"," ", OverlapID))%>%
  mutate(OverlapID= gsub("40","", OverlapID))%>%
  mutate(across(where(is.numeric), round, 2))%>%
  tibble::column_to_rownames("OverlapID")%>%
  as.data.frame()

#renames colnames
colnames(elp.olp_community)<- gsub(x =colnames(elp.olp_community), pattern = "_", replacement = " ")
colnames(elp.olp_community)<- gsub(x =colnames(elp.olp_community), pattern = "40", replacement = "")

#plot 
ggcorrplot::ggcorrplot(elp.olp_community, lab = T, outline.color = "white", lab_size = 3, tl.cex = 10)+
  scale_fill_gradient2(limit = c(0,1), low = "white", high = "grey50", mid = "grey80", midpoint = 0.5)+
  labs(fill="Overlap value")+
  theme(axis.text = element_text(face="italic", size = 10),
        legend.text = element_text(size=10),
        legend.title = element_text(size=10),
        plot.background = element_rect(colour = "white"))

Code 
#ggsave("matrix_overlap.png", path = "figures", dpi = 700, width = 8, height = 6)

3. Depth segregation

Cluster

  • input data for cluster = overlap matrix
  • depth distribution = from complete 2021 trawling data (not only individuals sampled for isotope)
Code 
# calculate median depth by species with trawling data set 2021
mean_depth_sp <- density_distribution%>%
  group_by(Nom_Scientifique)%>%
  mutate(median_depth=median(trawling_depth))%>%
  select(Nom_Scientifique, median_depth)%>%
  distinct()%>%
  ungroup()%>%
  arrange(Nom_Scientifique)%>%
  rename(species=Nom_Scientifique)

# calcul overlap ellipse 40%
test.elp_community <- rKIN::estEllipse(data= isotope_data_fish,  x="d13c", y="d15n", group="species", levels=40, smallSamp = TRUE)
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
[1] "data.frame"
Code 
elp.olp_community <- rKIN::calcOverlap(test.elp_community)

df <- elp.olp_community%>%
  tibble::column_to_rownames(var="OverlapID")%>%
  as.data.frame()

# Gap statistic
factoextra::fviz_nbclust(df, kmeans, nstart = 25,  method = "gap_stat", nboot = 50, verbose = FALSE)+
  labs(subtitle = "Gap statistic method")

Code 
res.km <- kmeans(scale(df),5, nstart = 25)

# Réduction de dimension en utilisant l'ACP
res.pca <- prcomp(df[, -5],  scale = TRUE)
# Coordonnées des individus
ind.coord <- as.data.frame(factoextra::get_pca_ind(res.pca)$coord)
# Ajouter les clusters obtenus à l'aide de l'algorithme k-means
ind.coord$cluster <- factor(res.km$cluster)
# Percentage of variance explained by dimensions
eigenvalue <- round(factoextra::get_eigenvalue(res.pca), 1)
variance.percent <- eigenvalue$variance.percent

# Add species names and median depth 
ind.coord_depth <- cbind(ind.coord, mean_depth_sp)
ind.coord_depth$median_depth <- as.factor(ind.coord_depth$median_depth)

Cluster and depth distribution

Code 
# Cluster plot ----
cluster_plot <- ggpubr::ggscatter(
  ind.coord_depth, x = "Dim.1", y = "Dim.2",
  color = "cluster", shape = "median_depth", 
  ellipse = TRUE, ellipse.type = "convex",
  palette = c("#007A75","#150578", "#B80050", "#845EC2", "#E77823"),
  size = 2.5, ggtheme = theme_minimal(),
  add.params = list(size=0.1),
  xlab = paste0("Dim 1 (", variance.percent[1], "% )" ),
  ylab = paste0("Dim 2 (", variance.percent[2], "% )" )) +
  scale_shape_manual(values = c(15, 8, 16, 2, 17, 18, 4))+
  labs(shape="Median depth (m)")+
  theme(plot.background = element_rect(color = "white"),
        legend.text = element_text(size=10.5),
        legend.title = element_text(size=10.5),
        axis.text = element_text(size=10.5),
        axis.title =element_text(size=10.5))+
  guides(col="none", fill="none")+
  ggrepel::geom_text_repel(label= ind.coord_depth$species, aes(col=cluster), size=4, box.padding = 0.5) +
  scale_color_manual(values = c("#007A75","#150578", "#B80050", "#845EC2", "#E77823"))

# Density plot----
# assigne each species to a cluster 
density_distribution_cluster <- density_distribution%>%
  mutate(cluster=case_when(Nom_Scientifique%in% c("Argyropelecus olfersii", "Lampanyctus crocodilus",
                                                  "Benthosema glaciale","Myctophum punctatum","Serrivomer beanii")~4,
                           Nom_Scientifique%in% c("Lampanyctus macdonaldi", "Maulisia argipalla",
                                                  "Searsia koefoedi")~5,
                           Nom_Scientifique%in% c("Cyclothone spp.","Notoscopelus bolini", "Notoscopelus kroyeri",
                                                  "Melanostigma atlanticum")~3,
                           Nom_Scientifique%in% c("Xenodermichthys copei", "Maurolicus muelleri")~1,
                           Nom_Scientifique%in% c("Arctozenus risso","Lestidiops sphyrenoides")~2))%>%
  group_by(Nom_Scientifique)%>%
  arrange(desc(trawling_depth))

# Order in function of median depth
density_distribution_cluster$Nom_Scientifique = with(density_distribution_cluster, reorder(Nom_Scientifique, cluster, min))  

density_plot <-ggplot(density_distribution_cluster,
                      aes(x = trawling_depth, y = Nom_Scientifique, group = Nom_Scientifique, 
                          col=factor(cluster), fill=factor(cluster)))+ 
  scale_fill_manual(values = c("#150578", "#E77823", "#007A75","#845EC2", "#B80050"))+
  scale_color_manual(values = c("#150578", "#E77823", "#007A75","#845EC2", "#B80050"))+
  ggridges::stat_density_ridges(quantile_lines = TRUE, quantiles = 0.5 , alpha=0.4, size=0.4)+
  theme_classic()+
  theme(axis.text.y = element_text(size=7))+
  scale_y_discrete(position = "right")+
  scale_x_reverse()+
  coord_flip()+
  ylab(label = "")+ xlab("Depth (m)")+
  theme(axis.text.y = element_text(size=10.5),
        axis.text.x = element_text(face="italic", size=10.5, angle=80, vjust = 0.5, hjust=0),
        axis.title.x = element_text(size=10.5),
        axis.title.y = element_text(size=10.5))+
  guides(fill="none", col="none", alpha="none")

ggpubr::ggarrange(cluster_plot, density_plot, ncol=1, labels = c("A", "B"), heights = c(1, 1.3))

Code 
#ggsave("cluster_density.png", path = "figures", dpi = 700, height = 8, width = 10)

4. Isotopic diversity index

Definitions

  • Isotopic divergence : tends to 1 when all points (or their weights) are located at the edge of the convex hull, i.e. when the oragnisms with the most extreme isotopic values dominate the food web

  • Isotopic dispersion: equal to 1 when most points (or their weights) are far from the center of gravity of the point group, i.e. when organisms tend to have contrasting isotopic values

  • Isotopic evenness: tends to 1 when all organisms are equally distributed in isotopic space

  • Isotopic uniqueness: tends to 1 when most organisms (or their weights) are isolated in isotopic space, i.e. when most organisms (or those with the highest abundance) are isolated in isotopic space, their isotopic values are very different from the rest of the organisms

Representativeness of sampling

  • Percentage of species biomass in each station
Code 
# Species % biomass and abundance by station 
species_abundance <- trawling_data_evhoe21%>%
  # selection of mesopelagic trawl 
  filter(Code_Station%in% c("Z0524", "Z0518", "Z0512", "Z0508", 
                            "Z0503", "Z0497", "Z0492"))%>%
  #deletion when only genus name and genus already sampled in isotopy + A. carbo
  filter(!Nom_Scientifique%in%c("Cyclothone braueri","Cyclothone microdon", 
                                "Myctophidae", "Aphanopus carbo"))%>%
  select(Code_Station, Nbr, Nom_Scientifique, Tot_V_HV)%>%
  group_by(Nom_Scientifique, Code_Station)%>%
  mutate(nb_ind= sum(Nbr))%>%
  select(-Nbr)%>%
  distinct()%>%
  group_by(Code_Station)%>%
  mutate(sum_biomass_station=sum(Tot_V_HV))%>%
  ungroup()%>%
  group_by(Code_Station, Nom_Scientifique)%>%
  mutate(pourcentage_biomass_sp= Tot_V_HV/sum_biomass_station*100)%>%
  select(Nom_Scientifique, pourcentage_biomass_sp, nb_ind, Code_Station)%>%
  mutate(across(where(is.numeric), round, 2))
  
htmltools::tagList(DT::datatable(species_abundance))

Calculation

  • Formatting of data and calculation of relative biomass within each station
Code 
# Load data 
isotope_data_fish <- isotope_data %>% filter (species != "Meganyctiphanes_norvegica")

# sourcing the R functions from 'si_div' R script
source("R/si_div.R")

# Format indiviudal_si
individuals_si <- isotope_data_fish %>%
  select(individual_code, station, d13c, d15n, species_code) %>%
  rename(indiv_ID=individual_code,
         d13C= d13c,
         d15N =d15n,
         Species_code= species_code)%>%
  arrange(Species_code)

# species_status_biomass ----
# with all species sampled (not only species sampled for isotope)

# load total trawling data evhoe 2021
trawling_data_evhoe21 <-  utils::read.csv(here::here("data", "trawling_data_evhoe_2021.csv"), sep = ";", header = T, dec = ".")

species_status_biomass <- trawling_data_evhoe21%>%
 # selection of mesopelagic trawl 
  filter(Code_Station%in% c("Z0524", "Z0518", "Z0512", "Z0508", 
                            "Z0503", "Z0497", "Z0492"))%>%
  #deletion when only genus name and genus already sampled in isotope + A. carbo
  filter(!Nom_Scientifique%in%c("Cyclothone braueri","Cyclothone microdon", "Myctophidae",
                                "Aphanopus carbo","Lampanyctus"))%>%
  select(Code_Station, Tot_V_HV, Nom_Scientifique, Code_Espece_Campagne)%>%
  distinct()%>%
  # sum species biomass by depth 
  group_by(Nom_Scientifique, Code_Station)%>%
  mutate(biomass_sp=sum(Tot_V_HV))%>%
  select(-Tot_V_HV)%>%
  distinct()%>%
  # sum of total biomass by station
  group_by(Code_Station)%>%
  mutate(biomass_tot=sum(biomass_sp))%>%
  # relative biomass of each species by station
  mutate(rel_biomass=biomass_sp/biomass_tot*100)%>%
  select(-c(biomass_sp, biomass_tot))%>%
  # selection of species sampled for isotopy in each depth 
  filter(Code_Station== "Z0492"& 
           Nom_Scientifique%in% c("Arctozenus risso", "Argyropelecus olfersii",
                                  "Lampanyctus crocodilus", "Myctophum punctatum",
                                  "Notoscopelus kroyeri", "Xenodermichthys copei")|
           Code_Station== "Z0497"& 
           Nom_Scientifique%in% c("Argyropelecus olfersii", "Lampanyctus crocodilus", 
                                  "Lampanyctus macdonaldi", "Myctophum punctatum",
                                  "Serrivomer beanii", "Xenodermichthys copei")|
           Code_Station== "Z0503"& 
           Nom_Scientifique%in% c("Arctozenus risso","Cyclothone","Lampanyctus crocodilus",
                                  "Serrivomer beanii", "Xenodermichthys copei")|
           Code_Station== "Z0508"& 
           Nom_Scientifique%in% c("Lestidiops sphyrenoides", "Maurolicus muelleri", 
                                   "Myctophum punctatum","Notoscopelus kroyeri",
                                  "Notoscopelus bolini")|
           Code_Station== "Z0512"& 
           Nom_Scientifique%in% c("Arctozenus risso","Argyropelecus olfersii",
                                  "Lampanyctus crocodilus","Notoscopelus kroyeri",
                                  "Xenodermichthys copei")|
           Code_Station== "Z0518"& 
           Nom_Scientifique%in% c("Argyropelecus olfersii","Lampanyctus crocodilus",
                                  "Benthosema glaciale", "Maulisia argipalla",
                                  "Searsia koefoedi", "Serrivomer beanii")|
           Code_Station== "Z0524"& 
           Nom_Scientifique%in% c("Argyropelecus olfersii","Lampanyctus crocodilus",
                                  "Melanostigma atlanticum",  "Serrivomer beanii",
                                  "Xenodermichthys copei"))%>%
  mutate(Species_code= tolower(Code_Espece_Campagne))%>%
  rename(Species_name= Nom_Scientifique,
         Status=Code_Station)%>%
  select(-Code_Espece_Campagne)%>%
  relocate(Status, .after=rel_biomass)%>%
  relocate(Species_code,.after = Species_name)%>%
  arrange(Species_name)

25m (Z0508)

Code 
# 25m Z0508 ----
individuals_si_25 <- individuals_si%>%
  filter(station =="Z0508")%>%
  select(-station)

status_biomass_25 <- species_status_biomass%>%
  filter(Status=="Z0508")

# computing mean Stable Isotope values for each species
# "group" column identical to species_code to fit with input format of function meanSI_group
# no "weight" input as number of indivuals sampled per species did not mirror actual species biomass
individuals_si_25<-data.frame(group=individuals_si_25[,"Species_code"], individuals_si_25)
mean_si_species_25<-meanSI_group(individuals_si_25)

# computing coefficent of variation within each species to assess intraspecific variability
cbind(CV_d13C=mean_si_species_25[,"sd_d13C"]/mean_si_species_25[,"d13C"], CV_d15N=mean_si_species_25[,"sd_d15N"]/mean_si_species_25[,"d15N"] )
               CV_d13C    CV_d15N
lesti-sph -0.012948117 0.03266115
maur-mue  -0.005224151 0.05249590
myct-pun  -0.009814716 0.04047466
noto-bol  -0.009178283 0.02755544
noto-kro  -0.012155251 0.02163464
Code 
# -> intraspecific variability is overall low (<20%)

# checking that species codes are the same in the two tables
#row.names(mean_si_species_25)==status_biomass_25[,"Species_code"] # OK

# building a single dataframe with all data for computing isotopic diversity indices
data_fish_25 <-data.frame(mean_si_species_25[,c("d13C","d15N", "sd_d13C","sd_d15N")], rel_biomass=status_biomass_25[,"rel_biomass"], Status=status_biomass_25[,"Status"], latin_name=status_biomass_25[,"Species_name"])

# scaling mean stable isotopes values using function "scale_rge01"
data_fish_scl_25<-scaleSI_range01(data_fish_25)

# computing isotopic diversity of the whole fish assemblage using scaled isotopic values and species relative biomass
ID_scl_ab_25<-IDiversity(cons=data_fish_scl_25, weight=data_fish_scl_25[,c("rel_biomass")], nm_plot="25m_Z0508")

# printing results
result <- as.data.frame(round(ID_scl_ab_25,3)) 
#knitr::kable(result)

Z0508_25m_d13C_d15N

370m (Z0492)

Code 
# 370m (Z0492) ----
individuals_si_370 <- individuals_si%>%
  filter(station=="Z0492")%>%
  select(-station)

status_biomass_370 <- species_status_biomass%>%
  filter(Status=="Z0492")

# computing mean Stable Isotope values for each species
# "group" column identical to species_code to fit with input format of function meanSI_group
# no "weight" input as number of indivuals sampled per species did not mirror actual species biomass
individuals_si_370<-data.frame(group=individuals_si_370[,"Species_code"], individuals_si_370)
mean_si_species_370<-meanSI_group(individuals_si_370)

# computing coefficent of variation within each species to assess intraspecific variability
cbind(CV_d13C=mean_si_species_370[,"sd_d13C"]/mean_si_species_370[,"d13C"], CV_d15N=mean_si_species_370[,"sd_d15N"]/mean_si_species_370[,"d15N"] )
              CV_d13C    CV_d15N
arct-ris -0.006981948 0.01886755
argy-olf -0.011380016 0.03779400
lamp-cro -0.004821553 0.04473307
myct-pun -0.019927081 0.04074265
noto-kro -0.012294978 0.02480797
xeno-cop -0.011716566 0.07257123
Code 
# -> intraspecific variability is overall low (<20%)

# checking that species codes are the same in the two tables
#row.names(mean_si_species_370)==status_biomass_370[,"Species_code"] # OK

# building a single dataframe with all data for computing isotopic diversity indices
data_fish_370 <-data.frame(mean_si_species_370[,c("d13C","d15N", "sd_d13C","sd_d15N")], rel_biomass=status_biomass_370[,"rel_biomass"], Status=status_biomass_370[,"Status"], latin_name=status_biomass_370[,"Species_name"])

# scaling mean stable isotopes values using function "scale_rge01"
data_fish_scl_370<-scaleSI_range01(data_fish_370)

# computing isotopic diversity of the whole fish assemblage using scaled isotopic values and species relative biomass
ID_scl_ab_370<-IDiversity(cons=data_fish_scl_370, weight=data_fish_scl_370[,c("rel_biomass")], nm_plot="370m_Z0492")

# printing results
result <- as.data.frame(round(ID_scl_ab_370,3)) 
#knitr::kable(result)

370m

555m (Z0512)

Code 
# 555m (Z0512) ----
individuals_si_555 <- individuals_si%>%
  filter(station=="Z0512")%>%
  select(-station)

status_biomass_555 <- species_status_biomass%>%
  filter(Status=="Z0512")

# computing mean Stable Isotope values for each species
# "group" column identical to species_code to fit with input format of function meanSI_group
# no "weight" input as number of indivuals sampled per species did not mirror actual species biomass
individuals_si_555<-data.frame(group=individuals_si_555[,"Species_code"], individuals_si_555)
mean_si_species_555<-meanSI_group(individuals_si_555)

# computing coefficent of variation within each species to assess intraspecific variability
cbind(CV_d13C=mean_si_species_555[,"sd_d13C"]/mean_si_species_555[,"d13C"], CV_d15N=mean_si_species_555[,"sd_d15N"]/mean_si_species_555[,"d15N"] )
              CV_d13C    CV_d15N
arct-ris -0.012564587 0.03643136
argy-olf -0.008743823 0.04560865
lamp-cro -0.017196928 0.05084072
noto-kro -0.014008937 0.01957619
xeno-cop -0.010054994 0.04599149
Code 
# -> intraspecific variability is overall low (<20%)

# checking that species codes are the same in the two tables
#row.names(mean_si_species_555)==status_biomass_555[,"Species_code"] # OK

# building a single dataframe with all data for computing isotopic diversity indices
data_fish_555 <-data.frame(mean_si_species_555[,c("d13C","d15N", "sd_d13C","sd_d15N")], rel_biomass=status_biomass_555[,"rel_biomass"], Status=status_biomass_555[,"Status"], latin_name=status_biomass_555[,"Species_name"])

# scaling mean stable isotopes values using function "scale_rge01"
data_fish_scl_555<-scaleSI_range01(data_fish_555)

# computing isotopic diversity of the whole fish assemblage using scaled isotopic values and species relative biomass
ID_scl_ab_555<-IDiversity(cons=data_fish_scl_555, weight=data_fish_scl_555[,c("rel_biomass")], nm_plot="555m_Z0512")

# printing results
result <- as.data.frame(round(ID_scl_ab_555,3)) 
#knitr::kable(result)

555m

715m (Z0503)

Code 
#715m Z0503 ----
individuals_si_715 <- individuals_si%>%
  filter(station=="Z0503")%>%
  select(-station)

status_biomass_715<- species_status_biomass%>%
  filter(Status=="Z0503")

# computing mean Stable Isotope values for each species
# "group" column identical to species_code to fit with input format of function meanSI_group
# no "weight" input as number of indivuals sampled per species did not mirror actual species biomass
individuals_si_715<-data.frame(group=individuals_si_715[,"Species_code"], individuals_si_715)
mean_si_species_715<-meanSI_group(individuals_si_715)

# computing coefficent of variation within each species to assess intraspecific variability
cbind(CV_d13C=mean_si_species_715[,"sd_d13C"]/mean_si_species_715[,"d13C"], CV_d15N=mean_si_species_715[,"sd_d15N"]/mean_si_species_715[,"d15N"] )
              CV_d13C    CV_d15N
arct-ris -0.013383792 0.02784671
cycl-otz -0.009264821 0.04946692
lamp-cro -0.013572275 0.03144289
serr-bea -0.008721937 0.06133097
xeno-cop -0.012017717 0.06196185
Code 
# -> intraspecific variability is overall low (<20%)

# checking that species codes are the same in the two tables
#row.names(mean_si_species_715)==status_biomass_715[,"Species_code"] # OK

# building a single dataframe with all data for computing isotopic diversity indices
data_fish_715 <-data.frame(mean_si_species_715[,c("d13C","d15N", "sd_d13C","sd_d15N")], rel_biomass=status_biomass_715[,"rel_biomass"], Status=status_biomass_715[,"Status"], latin_name=status_biomass_715[,"Species_name"])

# scaling mean stable isotopes values using function "scale_rge01"
data_fish_scl_715<-scaleSI_range01(data_fish_715)

# computing isotopic diversity of the whole fish assemblage using scaled isotopic values and species relative biomass
ID_scl_ab_715<-IDiversity(cons=data_fish_scl_715, weight=data_fish_scl_715[,c("rel_biomass")], nm_plot="715m_Z0503")

# printing results
result <- as.data.frame(round(ID_scl_ab_715,3)) 
#knitr::kable(result)

715m

1000m (Z0518)

Code 
#1000m Z0518 ----
individuals_si_1000 <- individuals_si%>%
  filter(station=="Z0518")%>%
  select(-station)

status_biomass_1000<- species_status_biomass%>%
  filter(Status=="Z0518")

# computing mean Stable Isotope values for each species
# "group" column identical to species_code to fit with input format of function meanSI_group
# no "weight" input as number of indivuals sampled per species did not mirror actual species biomass
individuals_si_1000<-data.frame(group=individuals_si_1000[,"Species_code"], individuals_si_1000)
mean_si_species_1000<-meanSI_group(individuals_si_1000)

# computing coefficent of variation within each species to assess intraspecific variability
cbind(CV_d13C=mean_si_species_1000[,"sd_d13C"]/mean_si_species_1000[,"d13C"], CV_d15N=mean_si_species_1000[,"sd_d15N"]/mean_si_species_1000[,"d15N"])
             CV_d13C    CV_d15N
argy-olf -0.01024511 0.03092443
bent-gla -0.01551299 0.06444602
lamp-cro -0.02274635 0.05835428
maul-arg -0.00960478 0.03191368
sear-koe -0.02285299 0.05388211
serr-bea -0.01480479 0.06708122
Code 
# -> intraspecific variability is overall low (<20%)

# checking that species codes are the same in the two tables
#row.names(mean_si_species_1000)==status_biomass_1000[,"Species_code"] # OK

# building a single dataframe with all data for computing isotopic diversity indices
data_fish_1000 <-data.frame(mean_si_species_1000[,c("d13C","d15N", "sd_d13C","sd_d15N")], rel_biomass=status_biomass_1000[,"rel_biomass"], Status=status_biomass_1000[,"Status"], latin_name=status_biomass_1000[,"Species_name"])

# scaling mean stable isotopes values using function "scale_rge01"
data_fish_scl_1000<-scaleSI_range01(data_fish_1000)

# computing isotopic diversity of the whole fish assemblage using scaled isotopic values and species relative biomass
ID_scl_ab_1000<-IDiversity(cons=data_fish_scl_1000, weight=data_fish_scl_1000[,c("rel_biomass")], nm_plot="1000m_Z0518")

# printing results
result <- as.data.frame(round(ID_scl_ab_1000,3)) 
#knitr::kable(result)

1000m

1335m (Z0497)

Code 
#1335 Z0497 ----
individuals_si_1335 <- individuals_si%>%
  filter(station=="Z0497")%>%
  select(-station)

status_biomass_1335 <- species_status_biomass%>%
  filter(Status=="Z0497")

# computing mean Stable Isotope values for each species
# "group" column identical to species_code to fit with input format of function meanSI_group
# no "weight" input as number of indivuals sampled per species did not mirror actual species biomass
individuals_si_1335<-data.frame(group=individuals_si_1335[,"Species_code"], individuals_si_1335)
mean_si_species_1335<-meanSI_group(individuals_si_1335)

# computing coefficent of variation within each species to assess intraspecific variability
cbind(CV_d13C=mean_si_species_1335[,"sd_d13C"]/mean_si_species_1335[,"d13C"], CV_d15N=mean_si_species_1335[,"sd_d15N"]/mean_si_species_1335[,"d15N"] )
              CV_d13C    CV_d15N
argy-olf -0.009889227 0.01043410
lamp-cro -0.018677334 0.05649516
lamp-mac -0.022058396 0.02734756
myct-pun -0.024242014 0.04136229
serr-bea -0.014132768 0.05883833
xeno-cop -0.012585074 0.05902294
Code 
# -> intraspecific variability is overall low (<20%)

# checking that species codes are the same in the two tables
#row.names(mean_si_species_1335)==status_biomass_1335[,"Species_code"] # OK

# building a single dataframe with all data for computing isotopic diversity indices
data_fish_1335 <-data.frame(mean_si_species_1335[,c("d13C","d15N", "sd_d13C","sd_d15N")], rel_biomass=status_biomass_1335[,"rel_biomass"], Status=status_biomass_1335[,"Status"], latin_name=status_biomass_1335[,"Species_name"])

# scaling mean stable isotopes values using function "scale_rge01"
data_fish_scl_1335<-scaleSI_range01(data_fish_1335)

# computing isotopic diversity of the whole fish assemblage using scaled isotopic values and species relative biomass
ID_scl_ab_1335<-IDiversity(cons=data_fish_scl_1335, weight=data_fish_scl_1335[,c("rel_biomass")], nm_plot="1335m_Z0497")

# printing results
result <- as.data.frame(round(ID_scl_ab_1335,3)) 
#knitr::kable(result)

bathypelagic

near_bottom (Z0524)

Code 
# 1010 Z0524 ----
individuals_si_nb <- individuals_si%>%
  filter(station=="Z0524")%>%
  select(-station)

status_biomass_nb <- species_status_biomass%>%
  filter(Status=="Z0524")

# computing mean Stable Isotope values for each species
# "group" column identical to species_code to fit with input format of function meanSI_group
# no "weight" input as number of indivuals sampled per species did not mirror actual species biomass
individuals_si_nb<-data.frame(group=individuals_si_nb[,"Species_code"], individuals_si_nb)
mean_si_species_nb<-meanSI_group(individuals_si_nb)

# computing coefficent of variation within each species to assess intraspecific variability
cbind(CV_d13C=mean_si_species_nb[,"sd_d13C"]/mean_si_species_nb[,"d13C"], CV_d15N=mean_si_species_nb[,"sd_d15N"]/mean_si_species_nb[,"d15N"] )
              CV_d13C    CV_d15N
argy-olf -0.016533711 0.03529939
lamp-cro -0.015886548 0.04471856
mela-atl -0.009854421 0.04078758
serr-bea -0.012700557 0.05654679
xeno-cop -0.009943149 0.04105076
Code 
# -> intraspecific variability is overall low (<20%)

# checking that species codes are the same in the two tables
#row.names(mean_si_species_nb)==status_biomass_nb[,"Species_code"] # OK

# building a single dataframe with all data for computing isotopic diversity indices
data_fish_nb <-data.frame(mean_si_species_nb[,c("d13C","d15N", "sd_d13C","sd_d15N")], rel_biomass=status_biomass_nb[,"rel_biomass"], Status=status_biomass_nb[,"Status"], latin_name=status_biomass_nb[,"Species_name"])

# scaling mean stable isotopes values using function "scale_rge01"
data_fish_scl_nb<-scaleSI_range01(data_fish_nb)

# computing isotopic diversity of the whole fish assemblage using scaled isotopic values and species relative biomass
ID_scl_ab_nb<-IDiversity(cons=data_fish_scl_nb, weight=data_fish_scl_nb[,c("rel_biomass")], nm_plot="near_bottom_Z0524")

# printing results
result <- as.data.frame(round(ID_scl_ab_nb,3)) 
#knitr::kable(result)

bathypelagic

Spider web

Code 
# load library
library(tidyr)

# crate a data frame with all index value
trophic_indices <- data.frame(
  group = c("25m", "370m", "555m", "715m", "1000m", "1335m", "1010m-bottom proximity"),
  IDiv = c(0.952, 0.918, 0.845, 0.928, 0.721, 0.932, 0.994),
  IDis = c(0.898, 0.620, 0.466, 0.621, 0.451, 0.755, 0.318),
  IEve= c(0.46, 0.713, 0.789, 0.644, 0.58, 0.633, 0.543), 
  IUni =c(0.566,0.884, 0.833, 0.599, 0.73, 0.713, 0.758))%>%
  mutate(across(group, factor, levels=c("25m", "370m", "555m", "715m", "1000m", "1335m", "1010m-bottom proximity")))

# plot 
ggradar::ggradar(trophic_indices,
                 values.radar = c("0", "0.5", "1"),
                 group.colours ="#004291", group.line.width = 1,
                 axis.label.size= 4, group.point.size = 3,
                 grid.label.size = 4)+
  facet_wrap(~group)+
  theme_minimal()+
  theme(strip.text.x = element_text(size=12, face = "bold"),
        plot.background = element_rect(colour = "white"))+
  guides(col="none")

Code 
#ggsave("spider_web.png", path = "figures", dpi = 700)

# spider web for each station 
# for (i in 1:7) {
#  print(ggradar::ggradar(trophic_indices[i,],
#                  values.radar = c("0", "0.5", "1"),
#                  group.colours ="#0c3d97",  group.line.width = 1.3,
#                  axis.label.size= 4, group.point.size = 5,
#                  grid.label.size = 4)+
#   facet_wrap(~group)+
#   theme_light()+
#   theme(strip.text.x = element_text(size=12, face = "bold"))+
#   guides(col="none"))
# }
htmltools::tagList(DT::datatable(trophic_indices))

5.additional analyses

Isotopic niches by sations

Code 
# colors selection 
nichecol <- c("#E4A33A", "#F67451", "#D664BE", "#3DA5D9", "#94b3ae",
              "#18206F", "#FD151B", "#049A8F", "#072AC8", "purple", 
              "#d193f7", "#d8c2ab", "#678FCB", "#A63A49", "#00547A", "#6B54A0")

# plot
ggplot(data = isotope_data_fish, 
                         aes(x = d13c, 
                             y = d15n)) + 
  facet_wrap(~factor(trawling_depth))+
  geom_point(aes(color = species), alpha=0.7, size=1) +
  scale_color_manual(values=nichecol)+  
  scale_fill_manual(values=nichecol)+
  scale_x_continuous(expression({delta}^13*C~'\u2030')) +
  scale_y_continuous(expression({delta}^15*N~'\u2030'))+
  stat_ellipse(aes(group = species, fill = species, color = species), 
               alpha = 0.2, level = p.ell,linewidth = 0.7, type = "norm", geom = "polygon")+
  theme_light()+
  theme(legend.text = element_text(size=10),
        legend.title = element_text(size=10),
        axis.title = element_text(size=10),
        axis.text = element_text(size=10))+
  labs(shape="Taxon", col= "Species", fill="Species")

Code 
#ggsave("niches_stations.png", path = "figures", dpi = 700, height = 8, width = 10)

Krill data

  • Total \(\delta\)15N variability = 0.76‰
  • Total \(\delta\)13C variability = 1.1‰
Code 
# Load data
isotope_data <-  utils::read.csv(here::here("data", "isotopic_data_2021.csv"), sep = ";", header = T,dec = ",")
isotope_data_krill <- isotope_data %>% filter (species == "Meganyctiphanes_norvegica")

first_plot <- ggplot(data = isotope_data_krill, aes(x = d13c, y = d15n)) + 
  geom_point(aes(color= factor(trawling_depth)), size = 3) +
  ylab(expression(paste(delta^{15}, "N (\u2030)"))) +
  xlab(expression(paste(delta^{13}, "C (\u2030)"))) + 
  theme(text = element_text(size=16)) + 
  theme_light()+
  paletteer::scale_color_paletteer_d("rcartocolor::Teal")

# Summarise By depth 
data_krill_sum <- isotope_data_krill %>% 
  group_by(trawling_depth) %>% 
  summarise(count = n(),
            mC = mean(d13c), 
            sdC = sd(d13c), 
            mN = mean(d15n), 
            sdN = sd(d15n))

second_plot <- first_plot +
  geom_errorbar(data = data_krill_sum,
                mapping = aes(x = mC, y = mN, ymin = mN - 1.96*sdN, ymax = mN + 1.96*sdN, col = factor(trawling_depth)), 
                width = 0, size=0.8) +
  geom_errorbarh(data = data_krill_sum, 
                 mapping = aes(x = mC, y = mN,xmin = mC - 1.96*sdC, xmax = mC + 1.96*sdC, col = factor(trawling_depth)),
                 height = 0, size=0.8) + 
  geom_point(data = data_krill_sum, 
             aes(x = mC, y = mN, fill = factor(trawling_depth),col = factor(trawling_depth), shape = factor(trawling_depth)),
             size = 5) +
  scale_shape_manual(values = c(21, 22, 23, 24, 25, 8, 7))+
  paletteer::scale_fill_paletteer_d("rcartocolor::Teal")+
  paletteer::scale_color_paletteer_d("rcartocolor::Teal")+
  labs(shape="Trawling depth (m)", col="Trawling depth (m)", fill= "Trawling depth (m)" )

print(second_plot)

Isotope variability

Code 
boxplot_d13c <- ggpubr::ggboxplot(isotope_data_krill , x = "trawling_depth", y = "d13c", 
          color = "trawling_depth",fill = "trawling_depth", alpha=0.3,
          ylab = "d13c", xlab = "trawling_depth")+
  ylab(expression(paste(delta^{13}, "C (\u2030)"))) +
  xlab("Trawling depth (m)")+
  paletteer::scale_fill_paletteer_d("rcartocolor::Teal")+
  paletteer::scale_color_paletteer_d("rcartocolor::Teal")+
  guides(fill="none", col="none" ,alpha="none")+
  theme_light()

boxplot_d15n <- ggpubr::ggboxplot(isotope_data_krill , x = "trawling_depth", y = "d15n", 
                  color = "trawling_depth", fill = "trawling_depth",
                  alpha=0.3) +
  ylab(expression(paste(delta^{15}, "N (\u2030)"))) +
  xlab("Trawling depth (m)")+
  paletteer::scale_fill_paletteer_d("rcartocolor::Teal")+
  paletteer::scale_color_paletteer_d("rcartocolor::Teal")+
  guides(fill="none", col="none" ,alpha="none")+
  theme_light()

ggpubr::ggarrange(boxplot_d13c, boxplot_d15n)

Tests

Code 
pairwise.wilcox.test(isotope_data_krill$d13c, isotope_data_krill$trawling_depth,
                 p.adjust.method = "BH")

    Pairwise comparisons using Wilcoxon rank sum exact test 

data:  isotope_data_krill$d13c and isotope_data_krill$trawling_depth 

     25    370   555   715   1000  1010 
370  0.359 -     -     -     -     -    
555  0.033 0.122 -     -     -     -    
715  0.042 0.125 0.497 -     -     -    
1000 0.033 0.033 0.303 0.179 -     -    
1010 0.033 0.033 0.541 0.125 0.883 -    
1335 0.107 0.433 0.833 1.000 0.454 0.433

P value adjustment method: BH
Code 
pairwise.wilcox.test(isotope_data_krill$d15n, isotope_data_krill$trawling_depth,
                 p.adjust.method = "BH")

    Pairwise comparisons using Wilcoxon rank sum exact test 

data:  isotope_data_krill$d15n and isotope_data_krill$trawling_depth 

     25    370   555   715   1000  1010 
370  0.025 -     -     -     -     -    
555  0.025 0.143 -     -     -     -    
715  0.025 0.243 0.030 -     -     -    
1000 0.025 0.274 0.585 0.025 -     -    
1010 0.025 0.120 0.629 0.025 0.467 -    
1335 0.025 0.306 0.037 0.750 0.025 0.025

P value adjustment method: BH